library(fixest) # to demean variables
library(multiwayvcov) # to cluster SEs
library(sandwich) # to get HC2 SEs
library(stargazer) # to generate tables
library(mvtnorm) # to simulate betas and vcov matrix

rm(list = ls())

# load data
load(file = 'chData.RData')

# Put each dataset in a different object
ch_reform <- ch_list[[1]]
ch_refTerm <- ch_list[[2]]
ch_pre_post <- ch_list[[3]]
ch_refTermMChanges <- ch_list[[4]]
ch_refTermParty <- ch_list[[5]]
ch_refTermMChangesParty <- ch_list[[6]]

# Exclude "Gonzalo Fuenzalida" and 'Javier Macaya' because their district changed after the reform
ch_reform <- subset(ch_reform, !(bill_author %in% c("Gonzalo Fuenzalida", 'Javier Macaya')))
ch_refTerm <- subset(ch_refTerm, !(bill_author %in% c("Gonzalo Fuenzalida", 'Javier Macaya')))
ch_pre_post <- subset(ch_pre_post, !(bill_author %in% c("Gonzalo Fuenzalida", 'Javier Macaya')))
ch_refTermMChanges <- subset(ch_refTermMChanges, !(bill_author %in% c("Gonzalo Fuenzalida", 'Javier Macaya')))
ch_refTermParty <- subset(ch_refTermParty, !(bill_author %in% c("Gonzalo Fuenzalida", 'Javier Macaya')))
ch_refTermMChangesParty <- subset(ch_refTermMChangesParty, !(bill_author %in% c("Gonzalo Fuenzalida", 'Javier Macaya')))

####### Models using M log ##########
ch_refTerm$lnmagWoman <- ch_refTerm$lnmag * ch_refTerm$female
ch_reform$lnmagWoman <- ch_reform$lnmag * ch_reform$female
ch_pre_post$lnmagWoman <- ch_pre_post$lnmag * ch_pre_post$female
ch_refTermMChanges$lnmagWoman <- ch_refTermMChanges$lnmag * ch_refTermMChanges$female

# M1: G*M C1,C2,C3 (all data)
data_m1 <- fixest::demean(women_pct + lnmag + female + lnmagWoman
	~ session, data = ch_refTerm)
data_m1 <- cbind(data_m1, bill_author = ch_refTerm$bill_author)
m1 <- lm(women_pct ~ lnmag + female + lnmagWoman - 1, data = data_m1)
vcov1 <- cluster.vcov(m1, cluster = data_m1$bill_author)
s1  <- sqrt(diag(vcov1))

# M2: G*M C1,C2,C3 (M changes at the middle of C2)
data_m2 <- fixest::demean(women_pct + lnmag + female + lnmagWoman
	~ session + reform, data = ch_refTermMChanges)
data_m2 <- cbind(data_m2, bill_author = ch_refTermMChanges$bill_author)
m2 <- lm(women_pct ~ lnmag + female + lnmagWoman - 1, data = data_m2)
vcov2 <- cluster.vcov(m2, cluster = data_m2$bill_author)
s2  <- sqrt(diag(vcov2))


# M3: G*M C1,C3
data_m3 <- fixest::demean(women_pct + lnmag + female + lnmagWoman
	~ session, data = subset(ch_refTerm, session %in% c('2010-2014', '2018-2022')))
data_m3 <- cbind(data_m3, bill_author = subset(ch_refTerm, session %in% c('2010-2014', '2018-2022'))$bill_author)
m3 <- lm(women_pct ~ lnmag + female + lnmagWoman - 1, data = data_m3)
vcov3 <- cluster.vcov(m3, cluster = data_m3$bill_author)
s3  <- sqrt(diag(vcov3))

# M4: C2
inSample <- names(table(ch_pre_post$bill_author)[table(ch_pre_post$bill_author) == 2])
ch_pre_post <- subset(ch_pre_post, bill_author %in% inSample)
data_m4 <- fixest::demean(women_pct + lnmag + female + lnmagWoman
	~ bill_author, data = ch_pre_post)
data_m4 <- cbind(data_m4, bill_author = ch_pre_post$bill_author)
m4 <- lm(women_pct ~ lnmag + lnmagWoman - 1, data = data_m4)
vcov4 <- cluster.vcov(m4, cluster = data_m4$bill_author)
s4  <- sqrt(diag(vcov4))

# M5: G*M C1,C2,C3 (only repeated legislators)
data_m5 <- fixest::demean(women_pct + lnmag + lnmagWoman
	~ bill_author + session, data = subset(ch_refTerm, allSessions == 1))
data_m5 <- cbind(data_m5, bill_author = subset(ch_refTerm, allSessions == 1)$bill_author)
m5 <- lm(women_pct ~ lnmag + lnmagWoman - 1, data = data_m5)
vcov5 <- cluster.vcov(m5, cluster = data_m5$bill_author)
s5  <- sqrt(diag(vcov5))

# M6: G*M C1, C2, C3 (M changes at the middle of C2) (only repeated legislators)
data_m6 <- fixest::demean(women_pct + lnmag + lnmagWoman
	~ session + reform + bill_author, data = subset(ch_refTermMChanges, allSessions == 1))
data_m6 <- cbind(data_m6, bill_author = subset(ch_refTermMChanges, allSessions == 1)$bill_author)
m6 <- lm(women_pct ~ lnmag + lnmagWoman - 1, data = data_m6)
vcov6 <- cluster.vcov(m6, cluster = data_m6$bill_author)
s6  <- sqrt(diag(vcov6))

# M7: G*M C1,C3 (only repeated legislators)
data_m7 <- fixest::demean(women_pct + lnmag + lnmagWoman
	~ bill_author + session, data = subset(ch_refTerm, session %in% c('2010-2014', '2018-2022') & l10_l18 == 1))
data_m7 <- cbind(data_m7, bill_author = subset(ch_refTerm, session %in% c('2010-2014', '2018-2022') & l10_l18 == 1)$bill_author)
m7 <- lm(women_pct ~ lnmag + lnmagWoman - 1, data = data_m7)
vcov7 <- cluster.vcov(m7, cluster = data_m7$bill_author)
s7  <- sqrt(diag(vcov7))

# M8: C3
data_m8 <- subset(ch_refTerm, session == '2018-2022')
m8 <- lm(women_pct ~ lnmag + female + lnmagWoman, data = data_m8)
vcov8 <- vcovHC(m8, 'HC2')
s8  <- sqrt(diag(vcov8))

# Table D.11
stargazer(m1, m2, m3, m4, m5, m6, m7, m8, 
	se = list(s1, s2, s3, s4, s5, s6, s7, s8), 
	covariate.labels = c('M (log)', 'Woman', 'M (log) x Woman'), 
	add.lines = list(c('FE by Legislative Term', 'Yes', 'Yes', 'Yes', 'No', 'Yes', 'Yes', 'Yes', 'No'),
		c('FE by Reform', 'No', 'Yes', 'No', 'No', 'No', 'Yes', 'No', 'No'),
		c('FE by Legislator', 'No', 'No', 'No', 'Yes', 'Yes', 'Yes', 'Yes', 'No')),
	dep.var.labels = "Bills' Portfolio on Women's Issues (%)",
	omit.stat = c("ser",'f'),
	star.cutoffs = c(0.10, 0.05, 0.01),
	no.space = TRUE,
	single.row = FALSE,	
	notes.append = FALSE,
	title = "Association Between Legislative Portfolio, Gender, and District Magnitude--2014-2021 Cámara de Diputadas y Diputados de Chile--Log Magnitude",
	label = 't:mag_log'
)

### Substantive Effect
## Scenarios: Women, M = 2
W2 <- data.frame(mag = log(2), female = 1, magAnalysisWoman = log(2) * 1)
## Scenarios: Women, M = 3 (post-reform)
W3 <- data.frame(mag = log(3), female = 1, magAnalysisWoman = log(3) * 1)
## Scenarios: Men, M = 2
M2 <- data.frame(mag = log(2), female = 0, magAnalysisWoman = log(2) * 0)
## Scenarios: Men, M = 3 (post-reform)
M3 <- data.frame(mag = log(3), female = 0, magAnalysisWoman = log(3) * 0)
# Women and Men Data
Wdata <- rbind(W2, W3)
Mdata <- rbind(M2, M3)

## Sims
set.seed(123)
sims1 <- rmvnorm(n = 1000, mean = coef(m1), sigma = vcov1)
sims2 <- rmvnorm(n = 1000, mean = coef(m2), sigma = vcov2)
sims3 <- rmvnorm(n = 1000, mean = coef(m3), sigma = vcov3)
sims4 <- rmvnorm(n = 1000, mean = coef(m4), sigma = vcov4)
sims5 <- rmvnorm(n = 1000, mean = coef(m5), sigma = vcov5)
sims6 <- rmvnorm(n = 1000, mean = coef(m6), sigma = vcov6)
sims7 <- rmvnorm(n = 1000, mean = coef(m7), sigma = vcov7)
sims8 <- rmvnorm(n = 1000, mean = coef(m8), sigma = vcov8)

## Calculate Predicted Values
# Women
pred1_W <- apply((sims1 %*% t(Wdata[2,])) - (sims1 %*% t(Wdata[1,])), 2, quantile, c(0.975, 0.5, 0.025))
pred2_W <- apply((sims2 %*% t(Wdata[2,])) - (sims2 %*% t(Wdata[1,])), 2, quantile, c(0.975, 0.5, 0.025))
pred3_W <- apply((sims3 %*% t(Wdata[2,])) - (sims3 %*% t(Wdata[1,])), 2, quantile, c(0.975, 0.5, 0.025))
pred4_W <- apply((sims4 %*% t(Wdata[2, -2])) - (sims4 %*% t(Wdata[1, -2])), 2, quantile, c(0.975, 0.5, 0.025))
pred5_W <- apply((sims5 %*% t(Wdata[2, -2])) - (sims5 %*% t(Wdata[1, -2])), 2, quantile, c(0.975, 0.5, 0.025))
pred6_W <- apply((sims6 %*% t(Wdata[2, -2])) - (sims6 %*% t(Wdata[1, -2])), 2, quantile, c(0.975, 0.5, 0.025))
pred7_W <- apply((sims7 %*% t(Wdata[2, -2])) - (sims7 %*% t(Wdata[1, -2])), 2, quantile, c(0.975, 0.5, 0.025))
pred8_W <- apply((sims8 %*% t(cbind(1, Wdata)[2,])) - (sims8 %*% t(cbind(1, Wdata)[1,])), 2, quantile, c(0.975, 0.5, 0.025))

# Men
pred1_M <- apply((sims1 %*% t(Mdata[2,])) - (sims1 %*% t(Mdata[1,])), 2, quantile, c(0.975, 0.5, 0.025))
pred2_M <- apply((sims2 %*% t(Mdata[2,])) - (sims2 %*% t(Mdata[1,])), 2, quantile, c(0.975, 0.5, 0.025))
pred3_M <- apply((sims3 %*% t(Mdata[2,])) - (sims3 %*% t(Mdata[1,])), 2, quantile, c(0.975, 0.5, 0.025))
pred4_M <- apply((sims4 %*% t(Mdata[2, -2])) - (sims4 %*% t(Mdata[1, -2])), 2, quantile, c(0.975, 0.5, 0.025))
pred5_M <- apply((sims5 %*% t(Mdata[2, -2])) - (sims5 %*% t(Mdata[1, -2])), 2, quantile, c(0.975, 0.5, 0.025))
pred6_M <- apply((sims6 %*% t(Mdata[2, -2])) - (sims6 %*% t(Mdata[1, -2])), 2, quantile, c(0.975, 0.5, 0.025))
pred7_M <- apply((sims7 %*% t(Mdata[2, -2])) - (sims7 %*% t(Mdata[1, -2])), 2, quantile, c(0.975, 0.5, 0.025))
pred8_M <- apply((sims8 %*% t(cbind(1, Mdata)[2,])) - (sims8 %*% t(cbind(1, Mdata)[1,])), 2, quantile, c(0.975, 0.5, 0.025))



# Building blocks of table D.12
bill_hat_women <- c(round(pred1_W[2,], 2), round(pred2_W[2,], 2), round(pred3_W[2,], 2),
	round(pred4_W[2,], 2), round(pred5_W[2,], 2), round(pred6_W[2,], 2),
	round(pred7_W[2,], 2), round(pred8_W[2,], 2))
ci_women <- c(paste0('(', round(pred1_W[3,], 2), ', ', round(pred1_W[1,], 2), ')'), 
	paste0('(', round(pred2_W[3,], 2), ', ', round(pred2_W[1,], 2), ')'), 
	paste0('(', round(pred3_W[3,], 2), ', ', round(pred3_W[1,], 2), ')'), 
	paste0('(', round(pred4_W[3,], 2), ', ', round(pred4_W[1,], 2), ')'), 
	paste0('(', round(pred5_W[3,], 2), ', ', round(pred5_W[1,], 2), ')'), 
	paste0('(', round(pred6_W[3,], 2), ', ', round(pred6_W[1,], 2), ')'), 
	paste0('(', round(pred7_W[3,], 2), ', ', round(pred7_W[1,], 2), ')'), 
	paste0('(', round(pred8_W[3,], 2), ', ', round(pred8_W[1,], 2), ')'))


bill_hat_men <- c(round(pred1_M[2,], 2), round(pred2_M[2,], 2), round(pred3_M[2,], 2),
	round(pred4_M[2,], 2), round(pred5_M[2,], 2), round(pred6_M[2,], 2),
	round(pred7_M[2,], 2), round(pred8_M[2,], 2))
ci_men <- c(paste0('(', round(pred1_M[3,], 2), ', ', round(pred1_M[1,], 2), ')'), 
	paste0('(', round(pred2_M[3,], 2), ', ', round(pred2_M[1,], 2), ')'), 
	paste0('(', round(pred3_M[3,], 2), ', ', round(pred3_M[1,], 2), ')'), 
	paste0('(', round(pred4_M[3,], 2), ', ', round(pred4_M[1,], 2), ')'), 
	paste0('(', round(pred5_M[3,], 2), ', ', round(pred5_M[1,], 2), ')'), 
	paste0('(', round(pred6_M[3,], 2), ', ', round(pred6_M[1,], 2), ')'), 
	paste0('(', round(pred7_M[3,], 2), ', ', round(pred7_M[1,], 2), ')'), 
	paste0('(', round(pred8_M[3,], 2), ', ', round(pred8_M[1,], 2), ')'))

tableC11 <- rbind("Women "= bill_hat_women, " " = ci_women, 
	"Men "= bill_hat_men, " " = ci_men)
colnames(tableC11) <- paste0('Model ', 1:8)
stargazer(tableC11,
	title = "Predicted percentage point change in the cosponsorship portfolio on women’s issues when M increases by one")

